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Abstract 



^ We study a family of distributions that arise in critical unitary random matrix 

ensembles. They are expressed as Fredholm determinants and describe the limiting 
^ | distribution of the largest eigenvalue when the dimension of the random matrices 

tends to infinity. The family contains the Tracy-Widom distribution and higher 
order analogues of it. We compute the distributions numerically by solving a 
Riemann-Hilbert problem numerically, plot the distributions, and discuss several 
properties that they appear to exhibit. 
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43 1 Introduction 

Consider the space of Hermitian nx n matrices with a probability measure of the form 

1 e -ntrV(M) dM ^ (L1) 
£ — Z n 

<N 

where dM is the Lebesgue measure defined by 

CO 

• n 

dM = n m h n dRe M ij dim M iji 

t— I 3=1 i<j 

and where the external field V is a real analytic function with sufficient growth at 
±oo. The quadratic case V(M) = \M 2 corresponds to the Gaussian Unitary Ensem- 
ble (GUE) and consists of matrices with independent Gaussian entries. The probability 
measure ( |1.1[ ) is invariant under conjugation with a unitary matrix and induces a prob- 
ability measure on the matrix eigenvalues given by 

Zn i<j j=l 

The large n limit of the average counting measure of the eigenvalues of a random 
matrix (hereafter referred to as the limiting mean eigenvalue distribution) exists and is 
characterized as the unique measure \xy which minimizes the logarithmic energy [7] 

MM) =l\ ^g — l —r dn{x)dM + [ V(x)dfx( X ), (1.3) 

JJ \x-y\ J 



1 



among all probability measures fi on M. The equilibrium measure \xy depends on the 
external field V and is supported on a finite union of intervals Uj =1 [aj, bj], with I, cij, bj 
depending on 7. It is absolutely continuous with a density tpy of the form [TU] 

I 

$vi?) = II \/ ( b J - x )( x - a j) h (x), x G \J j=1 [ aj ,bj], (1.4) 
i=i 

where h(x) is a real analytic function on R. Generically |14j . /i has no zeros on the 
support and in particular at the endpoints atj , bj , so that the equilibrium density has 
square root vanishing at the endpoints. However, there exist singular external fields V 
which are such that h(bj) = (or h(a,j) = 0). If h(bj) = 0, then necessarily an odd 
number of derivatives vanish: we have 

h'(bj) = ■■■ = h {2k - l) {bj) = 0, hW(bj) ^0, (1.5) 

for k G N. The value k = corresponds to the generic square-root behavior and k = 1 
to the first type of critical behavior where i^v{x) ~ c(bj — x) 5//2 as x /* bj. 

Example 1.1 The simplest critical case k = 1 occurs, e.g., for the critical quartic 
potential 

V(x) = ^ - ^x 3 + l -x 2 + |x. (1.6) 
The limiting mean eigenvalue density is then supported on [—2, 2] and given by [3] 

M*) = ^( x + 2) 1/2 ( 2 " *) 5/2 > x e t- 2 ' 2]- (1-7) 

Example 1.2 The general situation A; G N can be realized for instance with the fol- 
lowing family of polynomials, 

V k (x) - 2 • {2k + 2)l T (1 8) 

where (a) m = a • (a + 1) • • • (a + m — 1) is the Pochhammer symbol. The corresponding 
limiting mean eigenvalue densities are given by 

Mx) = ~ { fj, 2)l x^(l-x) 2fc+1 / 2 , xG [0,1]. (1.9) 



This can be verified by substituting (1.9) into the variational conditions that charac- 
terize the equilibrium measures in external field These conditions reduce to the 
equation 



Jo 



^ fc(s) ds = y fc / (x), for xG (0,1), 



x — s 

which is readily verified. 

The random matrix ensembles under consideration are known to have eigenvalues 
following a determinantal point process with a kernel of the form [7J 

e -2V(«) e -fV(y) 

K n (x,y) = (Pn(x)p n -l(y) -p n (y)Pn-l(x)), (1.10) 

x y n n 
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where pj is the degree j orthonormal polynomial with respect to the weight e~ nV ^ on 
M; Kj is its leading coefficient. Near a right endpoint b of the support for which k = 0, 
it was proved [5] that the large n limit of the re-scaled kernel is the Airy kernel: 

!™ -L K n (b+~^,b+ = K(°\u, v), (1.11) 



uniformly for u, t> > — L, L > 0, where c = an d -ft^ ) is the Airy kernel 

R{0) = Ai(K)Ai»-Ai(,)Ai'H (L12) 

u — u 

The limit of the probability that the largest eigenvalue is smaller than b + 2/3 is given 
by a Fredholm determinant: 

lim Prob (cn 2/3 (A n - 6) < s) = det(J - £T S (0) ), (1.13) 



where -fT^ denotes the Airy-kernel trace-class operator acting on L 2 (s, +00). The right 
hand side of the above equation is the Tracy-Widom distribution [19], which can be 
written in the form 

det(J - tf(°>) = exp (- £°°(y - s)q 2 (y)dy^ , (1.14) 

where go is the Hastings-McLeod solution go of the Painleve II equation 

q xx =xq + 2q 3 . (1.15) 

This solution [13] is characterized by the asymptotic conditions 

qo(x) ~ Ai (x), as x — > +00, (1-16) 

q (x) = J^- (1 + ^3 +0(x- 6 )Y asz^-oo. (1.17) 



If b is a right endpoint of the support for which k = 1, i.e. ( |1.5[ ) holds with = 1 
and 6j = b, the kernel has a different limit. It was showed in [BJ, see also 0], that 
there exists a limiting kernel K^' such that 

lim (1.18) 



n->oo cn 2 / 7 V cn 2 / 7 ' en 2 / 7 

uniformly for ii, t> in compact subsets of the real line. Since both sides of this equation 
tend to rapidly as u or v tends to +00, this result can be extended to a uniform 
statement for u, v > — L, L > 0, which implies 

lim Prob (cn 2/7 (A n - b) < s) = det(J - £T S (1) ), (1.19) 



where ifj is the trace class operator with kernel on (s,+oo). In more compli- 
cated double scaling limits, where V is not fixed but depends on n in a critical way, 
a limiting kernel K^\u, v ; to, ii) depending on two parameters was obtained. Those 
double scaling limits can describe phase transitions where the number of intervals in 
the support of the limiting mean eigenvalue distributions changes. If we consider an ex- 
ternal field V = Vt depending on a parameter t, it can happen that the support consists 
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of two intervals for t < 1, but of only one interval for t = 1. This happens for example 
when two intervals approach each other and simultaneously one of the intervals shrinks 
in the limit f\l. Such phase transitions are covered by the kernels K^\u, v;to,tx) 
with non-zero parameters tQ,t\. 

The kernels (u, v; to,ti) are related to the second member of the Painleve I hier- 
archy, but are most easily characterized in terms of a Riemann-Hilbert (RH) problem. 
We will give more details about this in the next section. 

Although no proofs have been given for k > 1, when considering in more detail the 
results and methods in [6] , one expects a result of the form 



m 2/(4 fc+ 3) ^ ( 6 + cn 2/lk + 3) > b + cn 2/(l + 3) ) = ^ *»' ■ ■ • > ^-l), (1-20) 



for general A; > 1, where the kernel now depends on 2/c parameters, and 

lim Prob (cn 2 /( 4fc+3 )(A n - 6) < a) = det(/ - #i fc) (*o, • • • ,*2fc-l))- (1.21) 



n— >oo 



The kernels if^ are related to the Painleve I hierarchy and will be characterized in 
the next section in terms of a RH problem. It was proved in [5] that the Fredholm 

determinant det(7 — Kg k \to, . . . , i2fc-i)) can be expressed explicitly in terms of a dis- 
tinguished solution to the equation of order 4fc + 2 in the second Painleve hierarchy, 
and in addition asymptotics for det(I — Kg k \to, . . . , i2fc-i)) as s — > ±oo were obtained. 
The asymptotics at +oo can be derived relatively easy from asymptotic properties of 
the kernel K^' and are given by 

4fc+3 

logdet(/-^ fc )(t ,...,t 2 fc-i)) = 0(e~ cs ^), ass^+oo. (1.22) 

The asymptotics as x — > — oo are more subtle and require a detailed analysis of the 
Fredholm determinants. In the simplest case to = • ■ • = £2^-1 = 0, they are given by 

3^2 



iogdet ( /-KW) = - * ON ^l^Lj *r +3 



4(4A: + 3)r(|) 2 r(2fe + 2) 21 
_ 2*+l lQg | s | +x (fc) + 0(\s\-^), ass -+-00, (1.23) 

o 

where T(x) is Euler's T-function. The constant has no explicit expression, except 
for k = 0, where it was proved in [l] that = 53 log 2 + £'(— 1), and £(s) is the 
Riemann zeta function. 

The goal of this paper is to set up a numerical scheme for computing the Fredholm 
determinants det(7 — Kg k \to, . . . , t^k— 1)), which will allow us to draw plots of the distri- 



butions and their densities, to verify formulas (1.22) and (1.23) numerically, to compute 
numerical values for the constants x \ an d to formulate a number of questions about 
the analytic properties of the distributions (monotonicity, inflection points), based on 
a closer inspection of the plots. In the next section, we define the kernels in a precise 
way using a RH problem. This RH characterization will also be used for the numerical 
analysis which we explain in more detail in Section |3l In Section [4] finally, we show 
plots of the distributions det(J — JQ ) and their densities for several values of k and 
the parameters to, ■ ■ ■ , i2fe-ij an d we will formulate a number of open problems. 
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Figure 1: The jump contour T and the jump matrices for 

2 Riemann— Hilbert characterization of the kernels 

The kernels K^> have the form 

^ ^ • • • ' ^ = -27T»(« - t,) ' (2 - 1} 

where the functions <&J- (to) = <&j\w\ to, ... , t 2 k-i) can be characterized in terms of 
a RH problem. 

RH problem for 

(a) $ = $( 2fc ) : C \ r -> C 2x2 is analytic, with 

r = uJ =1 rj-u{o}, ri = M + , r 3 = M _ , r 2 = e3S]R", r 4 = e^iR - , 

oriented as in Figure [T] 

(b) has continuous boundary values <£ + as £ approaches T \ {0} from the left, and 

from the right. They are related by the jump conditions 

<MC) = <MC)S;, forced (2.2) 

where 

Si=(l ]), (2-3) 



S2 = S,= (] J], (2.4) 



,0 1/ ' 

1 1) ' 
-1 

(c) <3> has the following behavior as £ — )• oo: 



= C~i CT3 iV f J + ha 3 C 1/2 + ©(C 1 )) e- 9 ^ 3 , (2.6) 
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2fc-l 



'1 x 
,0 -1 



where h = h(to, . . . , i2fc-i) is independent of £, 03 is the Pauli matrix 
iV is given by 

1 /1 1 



and 



»(«,...,M^(¥_ 2 g^, ( , 8) 

where the fractional powers are the principal branches analytic for £ £ C\(— 00, 0] 
and positive for £ > 0. 

(d) $ is bounded near 0. 

It was proved in [5] that this RH problem is uniquely solvable for any r eal values of 



(2k) (2k) 

to, . . . ,t2k-i- The functions $1 = &\ and $2 = ^2 appearing in (2.1) are the 



analytic extensions of the functions 3>n and $21 from the sector in between Ti and T2 
to the entire complex plane. Alternatively they can be characterized as fundamental 
solutions to the Lax pair associated to a special solution to the 2k-th member of the 
Painleve I hierarchy. We will not give details concerning this alternative description, 
since the RH characterization is more direct and more convenient for our purposes. 

Remark 2.1 The description in terms of differential equations in the PI hierarchy 
presents the possibility of computing these distributions using ODE solvers. However, 
similar to the Hastings-McLeod solution (see |18|). these solutions are inherently un- 
stable; hence, applying this approach in practice would require the use of high precision 
arithmetic, which is too computationally expensive to be practical. On the other hand, 
the representation in terms of a RH problem is numerically stable, and therefore is 
reliable. 



Not only the kernel can be described in terms of a RH problem, but also the 
logarithmic derivative of the Fredholm determinant can be expressed in terms of a RH 
problem, which shows similarities with the above one, but is nevertheless genuinely 
different. We have a formula of the form 

^ logdet(7 - #J*)(to, . . • ,t2*-i)) = 2^ (^ 1 (0^(C)) 21 l cv , (2-9) 

where X s is the unique solution to a RH problem, see jSJ Section 2]. 

This representation provides relative accuracy, whereas the representation as a Fred- 
holm determinant only provides absolute accuracy [2]. However, it requires solving a 
RH problem for each point of evaluation s and numerical indefinite integration to re- 
cover the distributions. Therefore, the expression in terms of a Fredholm determinant 
is more computationally efficient. 



3 Numerical study of the distributions 

We will compute the higher order Tracy- Widom distributions by calculating numer- 
ically, using the methodology of [16} [T5] . Consider the following canonical form for a 
RH problem: 
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Canonical form for RH problem for \& 

(a) \£ : C\r — > C 2x2 is analytic, where T is an oriented contour which is the closure of 
the set r = TiU - • - UT^ whose connected components can be Mobius-transformed 
to the unit interval Mi : Tj — > (—1, 1), with junction points r* = T \ T. 

(b) has continuous boundary values Vf^ as £ approaches T from the left, and 
from the right. For a given function G, they are related by the jump condition 

MC) = *-(OG(c). (3-i) 

(c) As £ — > oo, we have lim\I/(£) = /. 

(d) ^ is bounded near T* . 

Define the Cauchy transform 

and denote the limit from the left (right) for ( 6 T by Cjt (Cp ). We represent VP in 
terms of the Cauchy transform of an unknown function U defined on T: 

*(C) = / + C r t/(C). 



Plugging this into (3.1) we have the linear equation 



C^U -C r UG = G-I. (3.2) 

We solve this equation using a collocation method. 

We approximate U by U n for n = {n ri , . . . , vF l }, which is defined on each compo- 
nent Tj of the contour in terms of a mapped Chebyshev series: 

n r i — 1 

U n (x) = Uf'TjiMiix)), for x G I\ and i = 1, . . . , £, 
i=Q 

where G C 2x2 , and Tj is the j-th Chebyshev polynomial of the first kind. The 
convenience of this basis is that the Cauchy transforms Cp; [Tj o Mi] are known in closed 
form, in terms of hypergeometric functions which can be readily computed numerically 

D2|- 

For each £ G T*, let Ox, . . . , f2_L be the subset of components in T that have £ as an 
endpoint. In other words, Mj(£) = pj where pj = ±1 for i = 1, . . . , L. We say that U 
satisfies the zero sum condition if 

L 

J>i^(C) = o, 

i=l 

where U n ' denotes U restricted to Q{. The boundedness of if implies that U must 
satisfy the zero sum condition. 
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Define the mapped Chebyshev points of the first kind: 



/ 



Mr 1 



1 



\ 



COS 7T ( 1 



COS 



n i — 1 



V 



and the vector of unknown Chebyshev coefficients (in C 



2x2^ 



U 



WJeJ 

Then we can explicitly construct a matrix C~ such that 

/c r u n (^: 

crv = ; 

holds whenever U n satisfies the zero sum condition |15| . To define Cp C/ n (x ri ) at the 



C^U n (M^(±l))=hm i C^U n (Mr^ x )), 



endpoints, we use 

which exists when C7 n satisfies the zero sum condition. 
Thus we discretize (3.2) by 

L n U = (1 + C~)U - C-\JG n = G„ - / 

where G n = (G(x Fl ), . . . , G(x r *)) T and the multiplication by G n on the right is defined 
in the obvious way. 

The remarkable fact is that solving this linear system will generically imply that U n 
satisfies the zero sum condition if L n is nonsingular; if it does not, L n is necessarily not 
of full rank, and we can replace redundant rows with conditions imposing the zero sum 
condition |15j . Taking this possibly modified definition of L n , we have the following 
convergence result. 

Theorem 3.1 Ji5J / The L ro error of the numerical method is bounded by 

Cn 1 1 L n 1 1 oo \\\u-u n \\\, 

where C n grows logarithmically with max n ; U n is the polynomial which interpolates U 
at x Fl , . . . , x r * and 

III / III = WfWoo + max || (Mr 1 )'//]^. 



In practice, Hi" 1 1| appears to grow at most logarithmically with maxn whenever a 
solution to the RH problem exists. Therefore, if the solution U is smooth, the numerical 
method will converge spectrally as minn — > oo, with minn proportional to maxn. 

To apply the numerical method to the RH problem <£, we need to reduce it to 
canonical form. Define W{() = C~ aa/4: Ne- e ^ a,i , and we use the notation W± to 
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denote the analytic continuation of W above/below its branch cut along (—00, 0). We 
make the following transformation: 



(w(0 ICI>i 

1 I C| < 1 and ( lies between r 4 and T\ 

S± 1 S\ ICI < 1 and ( lies between Ti and r 2 ■ (3-3) 

S^SiS^ 1 |C| < 1 and ( lies between r 2 and T3 

/ |C| < 1 and £ lies between T3 and r 4 



Then $( 2fc ) satisfies the RH problem: 

RH problem for 

(a) f = ^( 2fe ) : C \ f -> C 2x2 is analytic, with 



r = ri u r 2 u r 4 u r 2 i u r 42 u r i4 u {i, , e«+3 }, 

fi = (l,oo), f 2 = -e4fcT3(oo, l), f 4 = -e4S+3(oo,l), 
f 21 = e i7r{1 ~«^' 0) ,f 42 = e^^^ifcTa^-SfeTs^fu = e i7r(0 

oriented as in Figure [2] 

(b) The jump conditions for ^ are given by 

MO = MCMO-W^tC), for C g fj, i = 1,2,4, 

tf+(C) = ^-(OS^W-Hc), for c e f 14 , 

*+(C) = v-iO^siw-Ho, for C g f 21 , 

MCHMOwr 1 ^), forcef 42 . 

(c) As £ — )• 00, we have lim^ r (^) = /. 

— 17T ITT 

(d) ^ is bounded near {1, e 4fc + 3 , e 4 *+ 3 }. 

With ^ in this form, we can readily compute it numerically, recover <I> by (3.3), and 



(k) 

thence evaluate the kernel of Kg numerically. This leaves one more task: computing 
the Fredholm determinant itself. We accomplish this using the framework of [2] , which 
also achieves spectral accuracy. 

4 Plots and open problems 

4.1 Local maxima of the densities 

In Figure [3j we plot the numerically computed distributions Ff.(s; 0, . . . , 0) for k = 
0, . . . , 5, where we write 

F k (s; t , ... , t 2 fc-i) = det(7 - A"i fc) (t , • • • , t 2k -i)). (4.1) 

In Figure [4] the corresponding densities are drawn. One observes that each of the 
densities has only one local maximum (i.e. the distributions have only one inflection 
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Figure 2: 



The jump contour T and the jump matrices for ^. 




-4-2 2 4 -6 -4 -2 2 4 



Figure 3: The distributions Fk for k = 0, 1, . . . , 5 with tj = 0. The slope steepens near 
— 1 when k increases. On the right, we also plot -Fqo (thick curve), as constructed in 
Section l4~4l 



point). The figures suggest that for any k £ N and for to = ... = t2fc-i = 0, the 
densities have only one local maximum. 

For general values of the parameters io)--->^2fc— l S M, the situation is different. 
We see in Figure [5j for k = 1, to = and varying negative ii, that the densities 
have two local maxima. From the random matrix point of view, this can be explained 
heuristically by the fact that the kernels K^'(u,v,Q,ti) for t\ < correspond to a 
double scaling limit which describes the transition from a random matrix model with 
a two-cut support (for the limiting mean eigenvalue distribution) to a one-cut support, 
where the parameter t\ regulates the speed of the transition. To be more precise, 
consider a random matrix ensemble with probability measure (1.1), where V = V n 
depends on n. If the dependence of V on n is fine-tuned in an appropriate way, it can 
happen that the equilibrium measure \iy n consists of two intervals for finite n, but of 
only one interval in the limit n — > oo. In order to obtain (u, v; 0, ii) as a scaling 
limit of the eigenvalue correlation kernel, both intervals in the support of [iy n should 
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-4-2 2 4 

Figure 4: The densities F' k (s) for k = 0, 1, . . . , 5 with i,- = 0. 
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0.5 1.0 1.5 2.0 2.5 3.0 3.5 

Figure 6: 1 — Ff.(s) for k = 1, . . . , 5 with tj = 0. 



approach each other and simultaneously one of the intervals should shrink, as n — > oo. 
If the n-dependence of V is chosen in an appropriate way, the limiting probability 
that a random matrix has an eigenvalue located in the shrinking interval lies strictly 
between and 1 (it actually increases when ti decreases). We believe that one local 
maximum of the densities in Figure [5] (the one most to the left) corresponds to the 
largest eigenvalue if no eigenvalues lie in the shrinking interval, and the second local 
maximum corresponds to the largest eigenvalue if this one lies in the shrinking interval. 
For k £ N, transitions can take place from at most k + 1 cuts to a one-cut regime, and 
for that reason we expect that for k E N, the density function has at most k + 1 local 
maxima, although we have no analytical evidence for this. 



4.2 Asymptotics as x — > +oo 

In Figure [6] we show the rate of convergence to one as s — > oo of -Ffc(s) for various values 
of k. For s < 1, we see that the distribution appears to approach a fixed distribution. 
For s > 1, the rate of convergence becomes increasingly rapid, matching the asymptotic 
formula (fl~22l). 



4.3 Asymptotics as x — > — oo 



We note that the constants i n (1-23) can be expressed as 
X {k) = lim ( logdet(/ - iff) - A?? 

s— >— OO V / 

cM 



lim ( logdet^-iT^) 



with 



A{V 



F(2k + 



3\2 
2) 



4(4fc + 3)r(|) 2 r(2£: + 2) 21 



d s \ogdet{I - K^)ds - Af^ 

|4fc+3 _ ^ + l log | s | 
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Error in approximating xq, scaled by [s| 

Error in approximating xo 




-250 -200 -150 -100 -50 -200 -150 -100 -50 



Figure 7: The absolute error in approximating xo, ~ logdet(7 — Kg ) + Ag\ as 
a function of s (left). The error multiplied by |s| 3 (right), showing faster convergence 
than predicted. 



For moderate M (we use M = —.5), we can reliably calculate logdet(J — K^) as 
before. For s < M, to reliably calculate d s logdet(I — Kg ), we use the RH problem for 
R used in Section 3.5]. This RH problem is in canonical form and d s logdet(/— Kg ) 
can be expressed in terms of its solution. We then expand d s logdet(I — Kg) — d s A s k ^ 
in piecewise Chebyshev polynomials, allowing for the efficient calculation of its integral. 

To verify the accuracy of the above approach, we need to estimate four errors, 
which we do using the following heuristics. We estimate the error in calculating vp( 2fe ) 
by ensuring that the smallest computed Chebyshev coefficient is below a given tolerance 
(10~ 12 ). The error in logdet(I — K^ff) is estimated by examining the Cauchy error as 
the number of quadrature points m in the Fredholm determinant routine increases. The 
error in <9 s logdet(I — K^) at each point of evaluation s is determined by examining 
the smallest computed Chebyshev coefficient of the numerical approximation to R. 
Finally, the accuracy of the piecewise Chebyshev approximation to d s log det(J — K^) 
is estimated by examining each piece's smallest Chebyshev coefficient. 

Using this approach, we estimate the first three 

•^(°) ~ —0.1365400105, (matches exact expression to 8 digits) 

-0.09614954 and 
X (2) ps -0.06145. 

Cancellation and other numerical issues cause the approach to be unreliable for larger 
k. 

The convergence to is verified in Figure [Tj One interesting thing to note is that 
the rate of convergence appears to be faster than predicted: numerical evidence suggest 
convergence like 0(|s| -3 ). A similar experiment for x^ suggests a convergence rate 
of 0(|s| -7 ). (The numerics for x^ are insuffici ently accurate to make a prediction. 



Therefore, we conjecture that the error term in (1.23) is in fact 0(\s\ ( 4fc + 3 )) 5 which is 
better than the theoretical error 0(\s\ 



4fc+3 , 
2 
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4.4 Large k limit 

For increasing k, one observes from Figure [3] that the slope of the distributions near 
— 1 gets steeper. At first sight, one may expect from Figure [3] that for large k, the 
distribution function tends to a step function, but a closer inspection reveals that this 
is not the case. Instead, we believe that there is a limit distribution supported on 
[—1,1] which is possibly discontinuous at —1 but continuous at 1. 

We present an asymptotic-numerical argument that this is indeed true. Consider 
the RH problem for ty( 2k \ Note that on fj, the jumps WSjW^ 1 — > I as k — > oo. Fur- 
thermore, inside the unit circle W(Q ~^ W^°°'(Q = (~ as ^N. Finally, T42 disappears. 
Thus, in a formal sense, we have the following RH problem: 



RH problem for #(°°) 
(a) * = ^r(oo) : C \ T ^ C 2x2 is analytic, with 

r = f 21 uf 14 u{±i}, f 21 = e^ 1 '°),f 14 



,i7r(Q,-l) 



(b) The jump conditions for \& are given by (for W = W^) 



* + (C) = *-(C)^ 1 5 1 w- 1 (C)> 



for C G T 14 , 
for C € f 21 . 



(c) As C — > 00, we have lim\l/(£) = /. 



This is not in canonical form: the jump matrices are not continuous at ±1, implying 
that the solution has singularities. We rectify this by using local parametrices to 
remove the jumps. Define 



h lo § (- 



■im - 1 

z — % 



2ni 



i. ° V z — % 



1 o 1 

} h 

1 1, 



|C- 1| < r, \z\ < 1, 

s 3 s 2 w-\ 

I C — 1 1 <r, \z\ > 1, 



and 



' e -2m/3 e 2i7r/3 N 
1 1 



' e -2i7r/3 e 2i7r/3 N 



-1 + i\ 



-1 + %\ 



0-3/6 / e - 2i7r /3 e 2i7r/3 N 
\ 1 1 



0-3/6 / e -2«7r/3 e 2«7r/3 



for |C + 1| < r, \z\ < 1, 
1 



1 



for |C + 1| < r, \z\ > 1, 



with the standard branch cuts, so that they lie on the half circle e^ ' t7T \ It is straight- 
forward to verify that P^ 1 ) have the same jumps as inside the disks |C -F 1| < r - 
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We now define for r sufficiently small 

no 



'#M(£)p(l) (£)-!, |C-l|<r, 
^)(C)p(-^)(C)-\ |C + l|<r, 
otherwise. 



Then, Y satisfies a RH problem in canonical form: 

RH problem for Y 

(a) Y : C \ A -> C 2x2 is analytic, with 

A = Ai U A 2 U r r (±l) U {±e ±ie }, 

Ax = e Ml-M) ! A 2 = e^H^- 1 ), r r (a) = {C : |C-«| =r}, 
where is given by r = \e l8 — 1|. 

(b) The jump conditions for Y are given by 

y+(C) = y-CO^r 1 ^^" 1 ^), for c e A l5 

y+(C) = y-to^iy- 1 ^), for c e a 2) 

y + (C) = y_(c)P (±1) (0, for C e r r (±i). 

(c) As £ — >• oo, we have iimy(£) = /. 

(d) y is bounded near {±e }. 

We can compute Y, and hence <I'( 00 ) numerically. We therefore define 



o |C|>1 

Sj 1 |C| < 1 and ImC < 
^S^S! |C| < 1 and ImC > 



which we use to compute the kernel of Ks°°^ , which is a trace-class operator now acting 
on L 2 (s, 1). (Again, we do not have a rigorous reason why the limiting operator acts 
only on L 2 (s, 1), not L 2 (s, oo). Instead, we justify this by the accuracy of the numerics.) 

It should be noted that the local parametrices P^ 1 ) are not close to the identity 
matrix on |C=p 1| = r, and therefore they would not be suitable parametrices to be used 
for a rigorous Deift/Zhou steepest descent analysis |12j applied to the RH problem 
for ^. However, this is not an issue here: numerically it is sufficient that the local 
parametrices satisfy the required jump conditions. 

While this construction has not been mathematically justified, it is perfectly usable 
in a numerical way. In fact, the resulting distribution matches the asymptotics for the 
finite k distributions, cf. Figure [3| providing strong evidence that, for — 1 < s < 1, 

lim det(/ - K^) = det(J - i^ (oo) ). 

We remark that the appears to be smooth near +1, hence Bornemann's numerical 
Fredholm determinant routine remains accurate for s > — 1. However, the singularity 
m at -1 causes the accuracy to break down as s approaches —1. Therefore, we 

cannot infer whether the distribution approaches zero smoothly, or if there is a jump. 
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